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Abstract: The interacting lattice gas model is used to simulate fluid flow 
through an open percolating porous medium with the fluid entering at the 
source-end and leaving from the opposite end. The shape of the steady-state 
concentration profile and therefore the gradient field depends on the porosity 
(p). The root mean square (rms) displacements of fluid and its constituents 
(tracers) show a drift power-law behavior, R oa t in the asymptotic regime 
(t oo). The flux current density (j) is found to scale with the porosity 
according to, j cx (Ap)^ with Ap = p — Pc and /? ~ 1.7. 

PACS numbers: 05.10.Ln; 05.40.-a; 0.5.60. Cd; 47.55.Mh; 83.85.Pt 
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1 Introduction 



Motion of fluid/gas constituents (tracers) determines the collective diffusion 
of fluid/gas concentration. In general, tracers movements are correlated and 
depend on the concentration of fluid and the porosity of host medium/matrix 
among other parameters such as temperature, interaction, pressure gradient, 
etc. Using a computer simulation model involving the interacting lattice gas, 
we consider the flow of fluid through a percolating porous medium [1] near 
its percolation threshold at a fixed temperature. It involves the motion of 
constituents, global transport of fluid, evolution of concentration profile in 
unsteady-flow, concentration gradient and flux rate in steady-state flow. A 
brief introduction and remark on some of these issues may help understanding 
our data. We simulate the flow through random walks (stochastic motion 
of each particle), but in contrast to many previous such diffusion studies we 
create a net flow by continuously injecting new particles at the bottom, thus 
forcing the already existing ones to move upwards since no two particles can 
occupy the same lattice site. The system thus is spatially inhomogeneous in 
vertical direction. 



1.1 Concentration diffusion 

Diffusion of fluid concentration (C) in a homogeneous space is described by 

§-o^v^c . (1) 

Evolution of the concentration profile, diffusion of its front in unsteady flow, 
and drift in steady-state flow regimes are well understood [2]. How does the 
concentration diffuse (i.e., fluid flow) in a porous medium [3] such as perco- 
lating system [4-6]? What type of motion do the fluid and its constituents 
exhibit in unsteady and steady-state flow regimes? How does the flux-rate 
depend on porosity? Addressing these questions becomes somewhat difficult 
with respect to solving the diffusion equation (1) particularly near perco- 
lation threshold [7-9] where the percolating pores are highly ramified; the 
boundary conditions involved with the pore space become prohibitively large 
to solve the diffusion equation numerically. We consider an interacting lattice 
gas to model the fiuid in order to address these questions in a percolating 
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porous matrix. The lattice gas consists of hard-core particles with a nearest 
neighbor interaction (see section 2). 

1.2 Tracer: diffusion, subdifFusion, drift 

A considerable progress has emerged in understanding the motion of a par- 
ticle executing its random walk in percolating system [4,5]. Variation of 
the root mean square displacement (R) with time (t) is one of the major 
quantities to characterize the type of motion, i.e., 

A-t"' + B -t"^ + ... (2) 

If Ui > 1/2, then Rc^ A-f^^ in the asymptotic time limit (t oo). With such 
a leading power- law variation, the motion is diffusive if i^i = 1/2 and drift if 
i/i — 1. In a percolating system at porosity (the fraction of pore sites) above 
threshold, p > Pc, we have fi — 1/2. At the percolation threshold [7-9], on 
the other hand, the random walk motion becomes anomalous diffusion with 
i/i ~ 0.2 in three dimensions and z/i ~ 0.3 in two dimensions [10]. 

In presence of a biased field [11-13], the motion of a particle can still be 

described by above equation (2) except at very high bias with diffusion (z/i = 
1/2) in short time and drift {v2 = 1) in the long time regime at p ^ pc with a 
crossover around t = {A/B)"^. The motion becomes very complex [11] as the 
porosity is reduced toward the threshold and the biased field competes with 
the barriers at the pore boundaries. The motion of a particle in such a porous 
medium has been extensively studied in a variety of biased fields for over a 
decade [11-17]. Many interesting findings have been reported particularly due 
to increased computing power. Some of these results include vanishing drift 
velocity above a characteristic or critical field [13,15], sub-diffusion and non- 
universal transport [11], log-periodic motion at high bias and large porosity 
[14]. Most of these studies deal with the motion of a single particle in a 
biased field. Instead, we consider the fiow of a fiuid of interacting particles 
through a percolating medium driven by concentration gradient (see below) . 
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1.3 Concentration gradient 



The mobile fluid constituents spread from high (toward the source end) to 
low particle concentration as the fluid flow from a source. The field caused 
by the concentration gradient drives the fluid constituents as described by 
the concentration diffusion equation (1) in a homogeneous space. Stream of 
particles emanating from the source execute their stochastic (random walk) 
motion. The instantaneous distribution of particles forms a special morphol- 
ogy. Particularly, the locus of the nearest neighbor particles on the moving 
fluid front, i.e., the profile of the concentration front, leads to the morphol- 
ogy of a percolating cluster at the percolation threshold. This is known as 
gradient percolation [18] since the dispersion and distribution of particles is 
caused by the concentration gradient. The shape of the concentration gradi- 
ent is well known and the concentration of particles on the front provides a 
good estimate of the percolation threshold. 

Evolution of the concentration profile of an interacting (nearest neighbor) 
lattice gas in a homogeneous space [19] seems consistent with the diffusion 
equation (1). In fact, the velocity of the front from eq. (1) can be used to cal- 
ibrate the time and length scale of the lattice gas simulations to understand 
the diffusion of specific system such as chlorine [19]. Obviously the motion 
of the front depends on the shape of the concentration profile, i.e., the con- 
centration gradient. How does the concentration profile change if the fluid 
moves through a percolating porous medium near its percolation threshold? 
How does the front speed depend on the porosity near the threshold? 

Let us consider a finite system with the one end (bottom) connected to the 
fluid source (particles) and the opposite end (top) open. The fiuid moves from 
bottom to top as driven by the concentration gradient. The fluid particles 
move from the source into the bottom and escape the system from the top. 
The concentration profile becomes stable as the system reaches the steady- 
state, i.e., when the in-flux (at the bottom) becomes equal to the out-fiux 
(from the top). The steady-state concentration profile depends on porosity 
and we investigate the changes in profile as we vary the porosity near the 
threshold. Since the shape of the profile, i.e. the concentration gradient, 
provides the driving field, the motion of the front and the fiux rate may 
depend on porosity as well. In this paper, we analyze some of these issues 
by a computer simulation model presented in next section (2) followed by 
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results and discussion (in section 3). The conclusion is provided in the last 
section 4. 



2 Model 

The host matrix is prepared on a simple-cubic L x L x L lattice. The per- 
colating porous medium of porosity p is generated by randomly distributing 
barriers, one at a site, on a fraction Pb = I — p sites. The porosity is kept 
above the percolation threshold p > Pc {— 0.312) so that a spanning cluster 
of connected pore extends from one of the sample to another. One end of this 
sample say the bottom (x = 0) is connected to a source of fluid represented 
by mobile particles. The opposite end (top) is open so that the fluid particles 
entering the lattice at the bottom can escape from the top. Each pore site in 
the bottom layer is occupied by the mobile fluid particle (with one particle 
at a site, the excluded volume effect). 

The fluid is modeled by an interacting lattice gas. In order to introduce 
interaction we assign an occupation variable {S) to each site. An empty pore 
site i is assigned Si — —1 while a site i with a fluid particle Si — 1 and with 
a barrier Si — 0. We use the following interaction energy, 

E = UY^S,S, (3) 

where the interaction strength U is set at unity in units of the Boltzmann 
constant {ks)- The summation is restricted to nearest neighbor sites in this 
simulation. The fluid-fluid repulsive and fluid-pore attractive interactions 
are thus considered. The fluid particles attempt to move to an empty neigh- 
bor in a randomly selected direction with a Boltzmann distribution as in 
the Metropohs algorithm at temperature T [20]. A periodic boundary con- 
dition is used along the transverse {y, z) directions and open condition along 
the longitudinal {x) direction: a fluid particle cannot move below the bot- 
tom plane while it can escape the system if it attempts to move above the 
top plane. As soon as a particle leaves the pore site at the bottom, it is 
occupied by another particle from the source. Thus a constant fluid concen- 
tration/density of unity is maintained in the bottom plane throughout the 
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simulation. An attempt to move each particles once is defined as one Monte 
Carlo step (MCS). 

As the simulation proceeds, particles move out into the medium from the 
bottom, fluid spreads, and concentration proflle evolves. The period during 
which none of the particles from the source reaches the top defines the un- 
steady ("short" time) regime. It is in this unsteady state regime that the 
lattice gas simulation seems to reproduce the results of continuum diffusion 
eq. (1) particularly the form of concentration proflle and its motion in homo- 
geneous space [19]. In the steady-state ("long" time) regime, the fluid-in-flux 
(at bottom) equals the fluid-out-flux (at the top), and the continuity equation 
for the conservation of the mass is satisfled, i.e., 

V-i = . (4) 

In the steady-state flow regime the concentration proflle becomes stable. One 
can evaluate the fluid current density (j) along the longitudinal direction. We 
investigate the transport behavior of fluid (the center of mass) and its con- 
stituent particles, the tracers and the flux response. Simulations are carried 
out for a long time with many independent runs to obtain a reliable estimate 
of physical quantities, particularly near percolation threshold. Further, we 
have used different lattice sizes to check for signiflcant flnite size effects which 
are not detected in these data. 



3 Results and Discussion 

Simulations are performed on different lattice sizes to look for flnite size 
effects with most data generated on 30^ and 50^ samples. The porosity is 
varied with p > 0.312. The temperature is constant T = 2 in units of 
Boltzmann constant. A^^ independent samples, K,. = 32 — 256, are used via 
parallelization with MPI calls with up to flve million time steps. We have 
analyzed the variation of the rms displacement of tracers and their center of 
mass with time, concentration proflle, and flux rate density as a function of 
porosity near percolation threshold. 



6 



3.1 RMS Displacements 



In section 1 we introduced the power-law dependence of the rms displacement 
(eq. 2) for a single particle executing its stochastic motion. The fluid consists 
of many particles (in our model) and the collective motion of the fluid (i.e., 
the center of mass) results from the motion of individual particles, the tracers. 
The rms displacement of tracer, R^, and that of the center of mass, Rc, are 
described by the power-law dependence (eq. 2). Since the behavior of Rc 
describes the motion of the fluid, the center of mass of the particles and 
fluid will be used synonymously as far as the fluid motion is concerned. As 
the fluid enters the system from the bottom, the concentration gradient field 
drive the fluid from the bottom. The longitudinal {x—) component of the 
rms displacement is much larger than the transverse {y, z) components. The 
power-law behavior of the total rms displacement (for both tracer and center 
of mass of the fluid) is dominated by the longitudinal component. Figure 
1 shows the variation of rms displacements for tracer with time at various 
porosities. The linear fit of data on a log-log scale in the short time regime 
(10^ — 10^ steps) suggest a power-law with exponent vi ~ 1/2. In the long 
time (asymptotic regime), the power-law exponent becomes higher V2 — 1 
leading to a drift-like motion. Such a crossover from diffusion to drift is 
clearly seen at relatively higher porosities {ps = 0.360,0.400). The crossover 
occurs at a characteristic time step (tc) which increases on decreasing the 
porosity toward the threshold. The rms displacement of the center of mass 
also shows similar crossover behavior (figure 2). 

The closer to the percolation threshold we are, the longer (tc) it takes to reach 
the asymptotic power-law behavior of the rms displacements for the tracer 
and fluid. Particularly at p = 0.312, the lowest porosity we simulated near 
percolation threshold {pc — 0.311608) [7], it was difficult to reach long time 
power-law behavior with the larger samples. Therefore, we performed our 
simulation on a smaller sample (30^) for 5 x 10® time steps. The variation of 
the longitudinal component of the rms displacement with time is presented 
in figure 3 which shows a clear crossover to drift in the long time regime. 
Further we note that the qualitative behavior for the variation of the rms 
displacement with time remain the same as that on the larger sample for 
both tracer and fluid. 

It is worth pointing out that there arc extensive computer simulations on 
biased diffusion with a variety of biased fields [11-17]. In most of these 



7 



studies the tracer's motion slows down at low porosity (closer to percolation 
threshold) and high bias as the bias competes with the barriers at the pore 
boundary. Thus there is no asymptotic drift behavior [15] of a biased ran- 
dom walk in a percolating medium near threshold at high biased field. In 
our driven system, the concentration gradient provides the bias. The gradi- 
ent field varies (see below) as we vary the porosity, but the asymptotic drift 
behavior of the rms displacement of the tracer and fluid persists no matter 
how close we are to the percolation threshold. This implies that the bias gen- 
erated by the concentration gradient, perhaps not too strong, is an effective 
method of transport through a porous medium. 



3.2 Concentration profile 



Initially, the concentration profile is a delta function at a; = 1, since there 
is no fluid particle in the sample except in the bottom plane {x = 1) where 
each pore site is occupied by a fluid particle. As the simulation proceeds the 
fluid spreads and the proflle changes. Figure 4 shows a typical evolution of 
the concentration profile. In the unsteady flow regime, i.e., when no fluid 
particle arrives at the top, the shape of the profile seems consistent with 
the concentration diffusion (eq. 1). In the long time (asymptotic) regime, 
the concentration profile becomes stable and the system reaches a steady- 
state fiow (see below). The steady-state profile depends on porosity. Figure 
5 shows the steady-state profile at various porosities above the percolation 
threshold. At low porosity (p = 0.312), we see a spatial variation in con- 
centration gradient. A relatively linear gradient (dc/dx) develops at higher 
porosities. The gradient field (concentration gradient) is delicately controlled 
by porosity in such a way that the asymptotic drift motion of tracer and fiuid 
occurs. The net transfer of fluid across the sample, i.e., the flow response, 
depends on porosity which is discussed in the following. 



3.3 Flow response 



Fluid enters the sample at the bottom {x — 1) and leaves from the top 
{x — L). There is a net fluid transfer [q) along x— direction. We evaluate 
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the current density (j) resulting from the rate of fluid (mass) transfer, 



1 dq 
L2 ■ It 



(5) 



In fact, we evaluate both, the input (bottom) and output (top) current den- 
sities (j) separately. Figure 6 shows the variation of these current (flux) 
densities at various porosities. As expected, the input flux density is higher 
and output flux density is lower initially, i.e. in the unsteady-flow regime. In 
steady-state, both must be equal, as seen in flgure 6. Thus we see that the 
flow has reached the steady-state at nearly all porosities except very close to 
the percolation threshold [ps = 0.312). Looking at the trend in data, it is 
easy to estimate the flux density by extrapolation. However, we have carried 
out o\ii simulation for a long time on a smaller sample to achieve the steady 
state at Ps — 0.312; flgure 7 shows the input and output flux density. 

It would be interesting to quantify how the steady-state flux density (j) 
depend on porosity. Figure 8 shows the variation of j with porosity (p^). We 
see a relatively good power-law dependence. 



where Ap = Ps —Pc and /3 ~ 1.7. Fig. 8 shows that this effective exponent (3 
increases with increasing L and thus asymptotically may be compatible with 
the exponent ~ 2 of random resistor networks [3-6] . 



4 Conclusion 

A computer simulation model is used to study the density proflle and fluid 
flow through an open porous medium near the percolation threshold. The 
fluid is driven by the concentration gradient which is evolved from maintain- 
ing a constant concentration, unity at the source (bottom) and zero at the 
top. RMS displacements of tracer and fluid are studied in detail at various 
porosities near the percolation threshold. The long time asymptotic power- 
law behavior of the rms displacement R oc t is found at all porosities above 
the threshold unlike the biased diffusion in previous simulations [11-15]. We 
believe that such an asymptotic drift is achieved due to the concentration 



j oc (Ap) 



(6) 



9 



gradient caused by a delicate competition between the stochastic motion of 
the fluid particles and the pore barriers. The shape of the density profile is 
sensitive to porosity. 

We have shown that the fluid flux through such a percolating porous medium 
reaches a steady-state above the percolation threshold. The flux density (j) 
depends on porosity and we characterize it by an empirical power-law relation 
(eq. 6). We plan to study the effects of parameters such as temperature, 
pressure, etc. in such a gradient driven system. 
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Figure 1: Root mean square (rms) displacement of the tracer with time 
steps on a log-log scale for various porosities, p = 0.312 — 0.4 on a 50^ sample 
with T = 2.0. Nr = 32 — 256 independent samples. 
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Figure 2: RMS displacement of the center of mass of the fluid particles with 
time steps. Statistics is the same as in figure 1. 
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Figure 3: RMS displacement of the tracer and the center of mass versus 
time steps at p = 0.312 on a 30^ sample with Nj. — 256. 
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Figure 5: Concentration profile in steady-state. Statistics is the same as in 
figure 1. 
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Figure 6: Flux rate density (j) versus time steps at various porosities. 
Statistics is the same as in figure 1. The top data is the flux-in at the 
bottom and bottom data is the flux-out from the top. 
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Figure 7: Flux rate density (j) versus time steps at p — 0.312 on a 30^ 
sample with the same statistics as in figure 3. 
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Figure 8: Log-log plot of j versus p — at T = 1 for L = 30 (x), 50 (+), 
60 (stars), and 70 (squares), with slopes 1.5, 1.6, 1.7, 1.75, respectively. We 
average over usually 8 lattices and 10^ time steps. This flux rate density j 
varies as the concentration gradient cal/L. 
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